publication: [IL-1β+ macrophages fuel pathogenic inflammation in pancreatic cancer] (https://www-nature-com.proxy.kib.ki.se/articles/s41586-023-06685-2#Sec1)

Published code on [GitHub:] (https://github.com/ostunilab/PDAC_Nature_2023/tree/main/scRNAseq/Mouse/Timecourse_KPC)

import libraries

library(Seurat)
## Loading required package: SeuratObject
## Loading required package: sp
## 'SeuratObject' was built under R 4.4.0 but the current version is
## 4.4.1; it is recomended that you reinstall 'SeuratObject' as the ABI
## for R may have changed
## 'SeuratObject' was built with package 'Matrix' 1.7.0 but the current
## version is 1.7.1; it is recomended that you reinstall 'SeuratObject' as
## the ABI for 'Matrix' may have changed
## 
## Attaching package: 'SeuratObject'
## The following objects are masked from 'package:base':
## 
##     intersect, t
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.5
## ✔ ggplot2   3.5.1     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.4     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(scDblFinder)
## Loading required package: SingleCellExperiment
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## 
## The following object is masked from 'package:dplyr':
## 
##     count
## 
## 
## Attaching package: 'MatrixGenerics'
## 
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## 
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## 
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## 
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## 
## The following object is masked from 'package:SeuratObject':
## 
##     intersect
## 
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## 
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, table,
##     tapply, union, unique, unsplit, which.max, which.min
## 
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## 
## The following objects are masked from 'package:lubridate':
## 
##     second, second<-
## 
## The following object is masked from 'package:tidyr':
## 
##     expand
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## 
## The following object is masked from 'package:utils':
## 
##     findMatches
## 
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## 
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## 
## The following object is masked from 'package:lubridate':
## 
##     %within%
## 
## The following object is masked from 'package:purrr':
## 
##     reduce
## 
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## 
## The following object is masked from 'package:sp':
## 
##     %over%
## 
## Loading required package: GenomeInfoDb
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## 
## Attaching package: 'Biobase'
## 
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## 
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
## 
## 
## Attaching package: 'SummarizedExperiment'
## 
## The following object is masked from 'package:Seurat':
## 
##     Assays
## 
## The following object is masked from 'package:SeuratObject':
## 
##     Assays
set.seed(123)

import data

dirs <-  list.dirs(path = './GSE217846_RAW/data/samples/', recursive = F, full.names = F)


for (x in dirs){
  
 cts <- ReadMtx(mtx= paste0('./GSE217846_RAW/data/samples/', x, '/matrix.mtx.gz'),
                             features = paste0('./GSE217846_RAW/data/samples/', x, '/features.tsv.gz'),
                             cells=paste0('./GSE217846_RAW/data/samples/', x, '/barcodes.tsv.gz'))
  
 assign(x, CreateSeuratObject(counts = cts, min.cells = 3, project = x))
 
}


# merge object

seurat_object <- merge(Healthy, y=c(Tumor_d10,Tumor_d20, Tumor_d30),
      add.cell.ids=c("Healthy","Tumor_d10", "Tumor_d20", "Tumor_d30"))

preprocessing

# create columns for metadata

# identity
unique(seurat_object@meta.data$orig.ident) # control identity
## [1] "Healthy"   "Tumor_d10" "Tumor_d20" "Tumor_d30"
# cell id
seurat_object@meta.data$cell_id <- rownames(seurat_object@meta.data)

# percent mt
seurat_object[["percent.mt"]]<-PercentageFeatureSet(seurat_object, pattern="^mt-")

plot1 <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

seurat_object <- subset(seurat_object, subset = percent.mt < 25 & nCount_RNA > 1000 & nFeature_RNA > 200)


head(seurat_object)
##                            orig.ident nCount_RNA nFeature_RNA
## Healthy_AAACCCAAGTAAACAC-1    Healthy       6472         1961
## Healthy_AAACCCACACTGCTTC-1    Healthy       6203         1950
## Healthy_AAACCCACAGCGACAA-1    Healthy       5495         1804
## Healthy_AAACCCATCGAAACAA-1    Healthy       3126         1095
## Healthy_AAACGAAAGACCTTTG-1    Healthy       1390          234
## Healthy_AAACGAAAGGAGTACC-1    Healthy       4835         1786
## Healthy_AAACGAAAGGAGTATT-1    Healthy       6673         2216
## Healthy_AAACGAAAGGCTAACG-1    Healthy       8308         2464
## Healthy_AAACGAACACTTGAAC-1    Healthy      10632         3141
## Healthy_AAACGAACAGGACAGT-1    Healthy       8371         2322
##                                               cell_id percent.mt
## Healthy_AAACCCAAGTAAACAC-1 Healthy_AAACCCAAGTAAACAC-1  14.292336
## Healthy_AAACCCACACTGCTTC-1 Healthy_AAACCCACACTGCTTC-1  13.557956
## Healthy_AAACCCACAGCGACAA-1 Healthy_AAACCCACAGCGACAA-1   7.770701
## Healthy_AAACCCATCGAAACAA-1 Healthy_AAACCCATCGAAACAA-1   8.925144
## Healthy_AAACGAAAGACCTTTG-1 Healthy_AAACGAAAGACCTTTG-1  10.503597
## Healthy_AAACGAAAGGAGTACC-1 Healthy_AAACGAAAGGAGTACC-1   7.321613
## Healthy_AAACGAAAGGAGTATT-1 Healthy_AAACGAAAGGAGTATT-1   6.069234
## Healthy_AAACGAAAGGCTAACG-1 Healthy_AAACGAAAGGCTAACG-1   7.029369
## Healthy_AAACGAACACTTGAAC-1 Healthy_AAACGAACACTTGAAC-1   6.593303
## Healthy_AAACGAACAGGACAGT-1 Healthy_AAACGAACAGGACAGT-1   9.090909
seurat_object
## An object of class Seurat 
## 22211 features across 54311 samples within 1 assay 
## Active assay: RNA (22211 features, 0 variable features)
##  4 layers present: counts.Healthy, counts.Tumor_d10, counts.Tumor_d20, counts.Tumor_d30
plot1 <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

# double check filtering
min(seurat_object@meta.data$nFeature_RNA)
## [1] 201
min(seurat_object@meta.data$nCount_RNA)
## [1] 1001
max(seurat_object@meta.data$percent.mt)
## [1] 24.99137

control for doublets

for (i in c('Healthy','Tumor_d10','Tumor_d20','Tumor_d30')){
    sub <- subset(seurat_object, subset = orig.ident == i)
    eval(parse(text=paste("sceDblF_",i," <- scDblFinder(GetAssayData(sub, layer = 'counts'), dbr = 0.07)",sep="")))
    eval(parse(text=paste("score.",i," <- sceDblF_",i,"@colData@listData[['scDblFinder.score']]",sep="")))
    eval(parse(text=paste("names(score.",i,") <- rownames(sceDblF_",i,"@colData)",sep="")))
}
## Assuming the input to be a matrix of counts or expected counts.
## Creating ~6189 artificial doublets...
## Dimensional reduction
## Evaluating kNN...
## Training model...
## iter=0, 1267 cells excluded from training.
## iter=1, 1263 cells excluded from training.
## iter=2, 1252 cells excluded from training.
## Threshold found:0.421
## 613 (7.9%) doublets called
## Assuming the input to be a matrix of counts or expected counts.
## Creating ~11676 artificial doublets...
## Dimensional reduction
## Evaluating kNN...
## Training model...
## iter=0, 2174 cells excluded from training.
## iter=1, 2442 cells excluded from training.
## iter=2, 2470 cells excluded from training.
## Threshold found:0.615
## 1779 (12.2%) doublets called
## Assuming the input to be a matrix of counts or expected counts.
## Creating ~12244 artificial doublets...
## Dimensional reduction
## Evaluating kNN...
## Training model...
## iter=0, 2233 cells excluded from training.
## iter=1, 2618 cells excluded from training.
## iter=2, 2544 cells excluded from training.
## Threshold found:0.532
## 1770 (11.6%) doublets called
## Assuming the input to be a matrix of counts or expected counts.
## Creating ~13341 artificial doublets...
## Dimensional reduction
## Evaluating kNN...
## Training model...
## iter=0, 2308 cells excluded from training.
## iter=1, 2665 cells excluded from training.
## iter=2, 2745 cells excluded from training.
## Threshold found:0.615
## 1995 (12%) doublets called
doublets.info <- rbind(sceDblF_Tumor_d10@colData,sceDblF_Tumor_d20@colData,sceDblF_Healthy@colData,sceDblF_Tumor_d30@colData)
seurat_object$is.doublet <- doublets.info$scDblFinder.class

seurat_object <- subset(seurat_object, subset = is.doublet == 'singlet')

seurat standard workflow

# default parameters
seurat_object<-NormalizeData(seurat_object)
## Normalizing layer: counts.Healthy
## Normalizing layer: counts.Tumor_d10
## Normalizing layer: counts.Tumor_d20
## Normalizing layer: counts.Tumor_d30
seurat_object<-FindVariableFeatures(seurat_object, selection.method = "vst", nfeatures=2000)
## Finding variable features for layer counts.Healthy
## Finding variable features for layer counts.Tumor_d10
## Finding variable features for layer counts.Tumor_d20
## Finding variable features for layer counts.Tumor_d30
all.genes<-rownames(seurat_object)
seurat_object<-ScaleData(seurat_object, features=all.genes)
## Centering and scaling data matrix
seurat_object <- RunPCA(seurat_object, features = VariableFeatures(object = seurat_object))
## PC_ 1 
## Positive:  Krt8, Krt18, S100a6, Clu, Lmna, Fbln2, Tm4sf1, Tm4sf4, Krt7, Krt19 
##     Sox4, Tspan8, Cavin1, Emp1, Cystm1, Onecut2, Nedd4, S100a16, 2200002D01Rik, Cdkn2a 
##     Cldn18, Anxa2, Scd2, Prkg2, Serpinb6a, Ccnd1, Tpm1, Pmepa1, Msln, Npnt 
## Negative:  Ctrb1, Try4, Prss2, Cela1, 2210010C04Rik, Cd74, Cela2a, Pnlip, Clps, Reg1 
##     Tyrobp, H2-Aa, Cpb1, Cela3b, Fcer1g, H2-Ab1, H2-Eb1, Try5, Ctrl, Cpa1 
##     Zg16, Rnase1, Pnliprp1, Sycn, Bcl2a1b, Il1b, Mpeg1, Ifi27l2a, Irf8, Gm2a 
## PC_ 2 
## Positive:  Tspan8, Krt18, Krt8, Cystm1, Krt7, Onecut2, Tm4sf4, Cldn18, Krt19, Lgals4 
##     Muc1, Cdh17, Lgals2, Cdkn2a, Msln, Ctse, Cdh1, Prxl2a, Prkg2, Apob 
##     Npnt, Anxa13, 2200002D01Rik, Cldn6, Krt20, Anxa10, Ccnd1, Dsp, Cldn2, Cpe 
## Negative:  Aebp1, Fstl1, Serping1, Bgn, Col3a1, Sparc, Col5a2, Loxl1, Dcn, Fbn1 
##     Col1a2, Meg3, Ccdc80, Col6a2, Serpinf1, Lhfp, Plod2, Adamts2, Slit3, C1s1 
##     Col1a1, Rcn3, Col6a3, Lox, Cpxm1, Gas1, Serpinh1, Prrx1, Cdh11, Pdgfra 
## PC_ 3 
## Positive:  Il1b, Lyz2, Cxcl2, Fcer1g, Ccl6, C1qa, Ptgs2, Mafb, Tgfbi, C1qc 
##     C1qb, Fcgr3, Csf1r, Tyrobp, Mpeg1, Thbs1, Apoe, Arg1, Dab2, Msr1 
##     Ccl9, Clec4n, Ctsd, Pla2g7, Cfp, Fn1, Ms4a6c, C3ar1, Trf, Csf2rb 
## Negative:  Cela2a, Cpa1, Cela3b, Pnlip, Ctrb1, Try4, Ctrl, Clps, Prss2, 2210010C04Rik 
##     Rnase1, Sycn, Cela1, Cpb1, Zg16, Reg1, Cel, Try5, Pnliprp1, Cpa2 
##     Pla2g1b, Uba52, Spink1, Tff2, Gimap4, Rps18, Klk1, Rpl21, Rps15a, Pnliprp2 
## PC_ 4 
## Positive:  Plvap, Emcn, Cdh5, Kdr, Egfl7, Mmrn2, Adgrf5, Esam, Pcdh17, Cyyr1 
##     Fabp4, Flt1, Tie1, Adgrl4, Myct1, Mcam, Ptprb, Prex2, Apold1, Gpihbp1 
##     Ramp2, Robo4, Epas1, Eng, Rasip1, Cd34, Fam167b, Grrp1, Arhgef15, Inhbb 
## Negative:  Dcn, Loxl1, Rpl32, Ccdc80, Serpinf1, Lox, Rpl21, Cpxm1, Col3a1, Col5a1 
##     Mfap5, Rps15a, Medag, Slit3, Pdgfra, Rps18, Rpl34, Col6a2, Lgi2, Cdh11 
##     Mmp2, Mmp23, Lgals1, Sfrp1, Fndc1, Rps2, Col5a2, Col6a1, Rcn3, Col1a2 
## PC_ 5 
## Positive:  Ambp, Cldn3, Dcdc2a, Tstd1, Epcam, Cldn7, Pigr, Ces1d, Onecut1, Prox1 
##     Ttr, Cp, Fxyd3, Ehf, Slc28a3, 5330417C22Rik, Smim22, Tmem229a, Atp1b1, Hpn 
##     Pkhd1, Slc5a1, Fgfr3, Wfdc15b, Nr5a2, Bche, Gsta3, Fmo2, Gstt1, Scara3 
## Negative:  Tmsb10, Top2a, F2r, Ube2c, Gimap4, Pclaf, Inpp4b, Thy1, Ptpn22, Tpx2 
##     Ccna2, Trbc2, Rrm2, Has2, Birc5, Smc4, Cd3g, Cd3e, Tuba1a, Ms4a4b 
##     Bcl2l15, Cenpe, Cenpa, S1pr1, Kif11, Cenpf, Hist1h2ap, Racgap1, Dynap, Hmgb2
ElbowPlot(seurat_object)

seurat_object <- RunUMAP(seurat_object, reduction='pca', dims = 1:20)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 10:35:50 UMAP embedding parameters a = 0.9922 b = 1.112
## 10:35:50 Read 48154 rows and found 20 numeric columns
## 10:35:50 Using Annoy for neighbor search, n_neighbors = 30
## 10:35:50 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 10:35:54 Writing NN index file to temp file /var/folders/8b/2_66jk7x0jl1qc4177kp8svm0000gq/T//RtmpdAYuqi/file474d3f45c7ff
## 10:35:54 Searching Annoy index using 1 thread, search_k = 3000
## 10:36:03 Annoy recall = 100%
## 10:36:07 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 10:36:08 Initializing from normalized Laplacian + noise (using RSpectra)
## 10:36:10 Commencing optimization for 200 epochs, with 2070706 positive edges
## 10:36:22 Optimization finished
seurat_object <- FindNeighbors(seurat_object, reduction = 'pca', dims = 1:20)
## Computing nearest neighbor graph
## Computing SNN
seurat_object <- FindClusters(seurat_object, resolution = c(0.4))
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 48154
## Number of edges: 1648610
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9476
## Number of communities: 23
## Elapsed time: 7 seconds
DimPlot(seurat_object, reduction='pca', label = TRUE)

DimPlot(seurat_object, label = TRUE)

Idents(seurat_object) <- 'RNA_snn_res.0.4'
p1 <- DimPlot(seurat_object, reduction = 'umap')
p2 <- DimPlot(seurat_object, reduction = 'umap',  group.by = c("orig.ident"))

p1+p2

integrate samples

# after preprocessing (standard seurat workflow)
seurat_object <- IntegrateLayers(object = seurat_object, method = CCAIntegration, assay = "RNA", orig.reduction = "pca", new.reduction = "integrated.cca",
                                 scale.layer = "scale.data",verbose = FALSE)

# follow up with dimensional reduction (again)
seurat_object <- RunUMAP(seurat_object, reduction='integrated.cca', dims = 1:20)
## 11:19:49 UMAP embedding parameters a = 0.9922 b = 1.112
## 11:19:49 Read 48154 rows and found 20 numeric columns
## 11:19:49 Using Annoy for neighbor search, n_neighbors = 30
## 11:19:49 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 11:19:51 Writing NN index file to temp file /var/folders/8b/2_66jk7x0jl1qc4177kp8svm0000gq/T//RtmpdAYuqi/file474d7cc535dd
## 11:19:52 Searching Annoy index using 1 thread, search_k = 3000
## 11:20:01 Annoy recall = 100%
## 11:20:02 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:20:03 Initializing from normalized Laplacian + noise (using RSpectra)
## 11:20:05 Commencing optimization for 200 epochs, with 2142582 positive edges
## 11:20:17 Optimization finished
seurat_object <- FindNeighbors(seurat_object, reduction = 'integrated.cca', dims = 1:20)
## Computing nearest neighbor graph
## Computing SNN
seurat_object <- FindClusters(seurat_object, resolution = c(0.6))
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 48154
## Number of edges: 1761701
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9197
## Number of communities: 23
## Elapsed time: 8 seconds
DimPlot(seurat_object, label = TRUE)

p1 <- DimPlot(seurat_object, reduction = 'umap', label = T)
p2 <- DimPlot(seurat_object, reduction = 'umap',  group.by = c("orig.ident"))

p1+p2

identify cells of interest

# fibroblast + CAF markers
VlnPlot(seurat_object, features = c("Col1a1", "Col1a2", "Col3a1",
                                    "Ngfr", "Pdgfra", "Acta2", "Cd74"))

FeaturePlot(seurat_object, features = c("Col1a1", "Col1a2", "Col3a1",
                                    "Ngfr", "Pdgfra", "Acta2", "Cd74", "Plvap"))

# macrophage markers
VlnPlot(seurat_object, features = c("C1qa", "Lyz2",
                                    "C1qb", "Mafb",
                                    "C1qc", "Apoe",
                                    "Csf1r", "Mrc1"))

FeaturePlot(seurat_object, features = c("C1qa", "Lyz2",
                                    "C1qb", "Mafb",
                                    "C1qc", "Apoe",
                                    "Csf1r", "Mrc1"))

# other immune cells

# monocytes
VlnPlot(seurat_object, features = c("Ccr2",  "Ly6c2"))

FeaturePlot(seurat_object, features = c("Ccr2",  "Ly6c2"))

# cDCs
VlnPlot(seurat_object, features = c("Flt3",  "Ccr7"))

FeaturePlot(seurat_object, features = c("Flt3",  "Ccr7"))

# neutrophils
VlnPlot(seurat_object, features = c("Cxcr2", "Csf3r"))

FeaturePlot(seurat_object, features = c("Cxcr2",  "Mpo", "Csf3r", "Elane"))

# pDCs
VlnPlot(seurat_object, features = c("Siglech",  "Ccr9"))

FeaturePlot(seurat_object, features = c("Siglech",  "Ccr9"))

# B cells
VlnPlot(seurat_object, features = c("Cd79b",  "Jchain", "Mzb1"))

FeaturePlot(seurat_object, features = c("Cd79b",  "Jchain", "Mzb1"))

# NK cells
VlnPlot(seurat_object, features = c("Gzma",  "Prf1"))

FeaturePlot(seurat_object, features = c("Gzma",  "Prf1"))

# T cells
VlnPlot(seurat_object, features = c("Lef1",  "Il2rb", "Cd3e", "Gzmb", "Cd8a", "Xcl1"))

FeaturePlot(seurat_object, features = c("Lef1",  "Il2rb", "Cd3e", "Gzmb", "Cd8a", "Xcl1"))

# tumor markers
VlnPlot(seurat_object, features = c("Epcam",  "Krt8", "Krt7", "Krt19",  "Top2a", "Gata6", "Hmga2"))

FeaturePlot(seurat_object, features = c("Epcam",  "Krt8", "Krt7", "Krt19",  "Top2a", "Hmga2", "Gata6"))

# acinar, ductal and ADM
VlnPlot(seurat_object, features = c("Krt19", "Cftr",  "Cpb1", "Cpa1", "Cela2a", "Sox9", "Spp1", "Reg3a", "Crp"))

FeaturePlot(seurat_object, features = c("Krt19", "Cftr",  "Cpb1", "Cpa1", "Cela2a", "Sox9", "Spp1", "Reg3a", "Crp"))

# endothelial cells
VlnPlot(seurat_object, features = c("Plvap",  "Cdh5"))

FeaturePlot(seurat_object, features = c("Plvap",  "Cdh5"))

subset to fibroblasts, tumor, acinar and epithelial cells

# rename identified clusters
seurat_object <- RenameIdents(object = seurat_object, `10` = "Fibroblasts_1")
seurat_object <- RenameIdents(object = seurat_object, `14` = "Fibroblasts_2")
seurat_object <- RenameIdents(object = seurat_object, `20` = "Fibroblasts_3")

seurat_object <- RenameIdents(object = seurat_object, `1` = "Macrophages_1") # includes monocytes
seurat_object <- RenameIdents(object = seurat_object, `2` = "Macrophages_2")
seurat_object <- RenameIdents(object = seurat_object, `16` = "Macrophages_3")

seurat_object <- RenameIdents(object = seurat_object, `15` = "Ductal/ADM cells")

seurat_object <- RenameIdents(object = seurat_object, `13` = "Acinar cells")

seurat_object <- RenameIdents(object = seurat_object, `0` = "Tumor cells_1")
seurat_object <- RenameIdents(object = seurat_object, `4` = "Tumor cells_2")
seurat_object <- RenameIdents(object = seurat_object, `8` = "Tumor cells_3")
seurat_object <- RenameIdents(object = seurat_object, `9` = "Tumor cells_4")

seurat_object <- RenameIdents(object = seurat_object, `11` = "cDCs")
seurat_object <- RenameIdents(object = seurat_object, `21` = "pDCs")

seurat_object <- RenameIdents(object = seurat_object, `5` = "Neutrophils_1")
seurat_object <- RenameIdents(object = seurat_object, `18` = "Neutrophils_2")
seurat_object <- RenameIdents(object = seurat_object, `22` = "Neutrophils_3")

seurat_object <- RenameIdents(object = seurat_object, `3` = "Bcells_1")
seurat_object <- RenameIdents(object = seurat_object, `19` = "Bcells_2")

seurat_object <- RenameIdents(object = seurat_object, `6` = "NKcells")

seurat_object <- RenameIdents(object = seurat_object, `7` = "Tcells")

seurat_object <- RenameIdents(object = seurat_object, `12` = "Endothel")

DimPlot(seurat_object, reduction = "umap", label = T) 

# save for easier import
saveRDS(seurat_object, file = "./GSE217846_RAW/data/timecourseKPC_preprocessed.rds")
sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: aarch64-apple-darwin20
## Running under: macOS 15.4.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Europe/Stockholm
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] scDblFinder_1.18.0          SingleCellExperiment_1.26.0
##  [3] SummarizedExperiment_1.34.0 Biobase_2.64.0             
##  [5] GenomicRanges_1.56.2        GenomeInfoDb_1.40.1        
##  [7] IRanges_2.38.1              S4Vectors_0.42.1           
##  [9] BiocGenerics_0.50.0         MatrixGenerics_1.16.0      
## [11] matrixStats_1.5.0           lubridate_1.9.3            
## [13] forcats_1.0.0               stringr_1.5.1              
## [15] purrr_1.0.4                 readr_2.1.5                
## [17] tidyr_1.3.1                 tibble_3.2.1               
## [19] ggplot2_3.5.1               tidyverse_2.0.0            
## [21] dplyr_1.1.4                 Seurat_5.2.1               
## [23] SeuratObject_5.0.2          sp_2.2-0                   
## 
## loaded via a namespace (and not attached):
##   [1] RcppAnnoy_0.0.22          splines_4.4.1            
##   [3] later_1.4.1               BiocIO_1.14.0            
##   [5] bitops_1.0-9              polyclip_1.10-7          
##   [7] XML_3.99-0.17             fastDummies_1.7.5        
##   [9] lifecycle_1.0.4           edgeR_4.2.2              
##  [11] globals_0.16.3            lattice_0.22-6           
##  [13] MASS_7.3-61               magrittr_2.0.3           
##  [15] limma_3.60.6              plotly_4.10.4            
##  [17] sass_0.4.9                rmarkdown_2.29           
##  [19] jquerylib_0.1.4           yaml_2.3.10              
##  [21] metapod_1.12.0            httpuv_1.6.15            
##  [23] sctransform_0.4.1         spam_2.11-1              
##  [25] spatstat.sparse_3.1-0     reticulate_1.40.0        
##  [27] cowplot_1.1.3             pbapply_1.7-2            
##  [29] RColorBrewer_1.1-3        abind_1.4-8              
##  [31] zlibbioc_1.50.0           Rtsne_0.17               
##  [33] RCurl_1.98-1.16           GenomeInfoDbData_1.2.12  
##  [35] ggrepel_0.9.6             irlba_2.3.5.1            
##  [37] listenv_0.9.1             spatstat.utils_3.1-2     
##  [39] goftest_1.2-3             RSpectra_0.16-2          
##  [41] dqrng_0.4.1               spatstat.random_3.3-2    
##  [43] fitdistrplus_1.2-2        parallelly_1.42.0        
##  [45] DelayedMatrixStats_1.26.0 codetools_0.2-20         
##  [47] DelayedArray_0.30.1       scuttle_1.14.0           
##  [49] tidyselect_1.2.1          UCSC.utils_1.0.0         
##  [51] farver_2.1.2              viridis_0.6.5            
##  [53] ScaledMatrix_1.12.0       spatstat.explore_3.3-4   
##  [55] GenomicAlignments_1.40.0  jsonlite_1.9.0           
##  [57] BiocNeighbors_1.22.0      progressr_0.15.1         
##  [59] scater_1.32.1             ggridges_0.5.6           
##  [61] survival_3.7-0            tools_4.4.1              
##  [63] ica_1.0-3                 Rcpp_1.0.14              
##  [65] glue_1.8.0                gridExtra_2.3            
##  [67] SparseArray_1.4.8         xfun_0.51                
##  [69] withr_3.0.2               fastmap_1.2.0            
##  [71] bluster_1.14.0            digest_0.6.37            
##  [73] rsvd_1.0.5                timechange_0.3.0         
##  [75] R6_2.6.1                  mime_0.12                
##  [77] colorspace_2.1-1          scattermore_1.2          
##  [79] tensor_1.5                spatstat.data_3.1-4      
##  [81] generics_0.1.3            data.table_1.16.4        
##  [83] rtracklayer_1.64.0        httr_1.4.7               
##  [85] htmlwidgets_1.6.4         S4Arrays_1.4.1           
##  [87] uwot_0.2.2                pkgconfig_2.0.3          
##  [89] gtable_0.3.6              lmtest_0.9-40            
##  [91] XVector_0.44.0            htmltools_0.5.8.1        
##  [93] dotCall64_1.2             scales_1.3.0             
##  [95] png_0.1-8                 spatstat.univar_3.1-1    
##  [97] scran_1.32.0              knitr_1.49               
##  [99] rstudioapi_0.17.1         tzdb_0.4.0               
## [101] reshape2_1.4.4            rjson_0.2.23             
## [103] nlme_3.1-166              curl_6.2.1               
## [105] cachem_1.1.0              zoo_1.8-12               
## [107] KernSmooth_2.23-24        vipor_0.4.7              
## [109] parallel_4.4.1            miniUI_0.1.1.1           
## [111] ggrastr_1.0.2             restfulr_0.0.15          
## [113] pillar_1.10.1             grid_4.4.1               
## [115] vctrs_0.6.5               RANN_2.6.2               
## [117] promises_1.3.2            BiocSingular_1.20.0      
## [119] beachmat_2.20.0           xtable_1.8-4             
## [121] cluster_2.1.6             beeswarm_0.4.0           
## [123] evaluate_1.0.3            locfit_1.5-9.11          
## [125] cli_3.6.4                 compiler_4.4.1           
## [127] Rsamtools_2.20.0          rlang_1.1.5              
## [129] crayon_1.5.3              future.apply_1.11.3      
## [131] labeling_0.4.3            ggbeeswarm_0.7.2         
## [133] plyr_1.8.9                stringi_1.8.4            
## [135] viridisLite_0.4.2         deldir_2.0-4             
## [137] BiocParallel_1.38.0       munsell_0.5.1            
## [139] Biostrings_2.72.1         lazyeval_0.2.2           
## [141] spatstat.geom_3.3-5       Matrix_1.7-1             
## [143] RcppHNSW_0.6.0            hms_1.1.3                
## [145] patchwork_1.3.0           sparseMatrixStats_1.16.0 
## [147] future_1.34.0             statmod_1.5.0            
## [149] shiny_1.10.0              ROCR_1.0-11              
## [151] igraph_2.1.4              bslib_0.9.0              
## [153] xgboost_1.7.8.1